function [ExoTimeStates,RegionsForExoStates,pproc_boot] = GenExoTimeStatesBoot(DS,opt,pproc)
% Create the time varying exogenous state variables indices
% pproc is the key 

% Load price data:

load([opt.DataFolder opt.PricesData])

% Set ARIMA model for sugar price:
MdlPS = arima(1,0,0);

% Get sugar price used:
sug = double(prices(:,opt.SugarPrice));

T = sum(~isnan(sug));
bootsugar(1) = pproc.sugar.Cte/(1-pproc.sugar.AR);

% Generate time series of sugar price based on estimated process:
for t=2:300
    bootsugar(t) = pproc.sugar.Cte + pproc.sugar.AR*bootsugar(t-1) + randn(1,1)*pproc.sugar.Sigma;
end
bootsugar = bootsugar(end-T+1:end);

% Estimate AR(1) for sugar:

%disp(' Estimate process for sugar prices: ')
MdlPS = estimate(MdlPS,bootsugar','Display','off');

pproc_boot.sugar.Cte    = MdlPS.Constant; % This is the the constant.
pproc_boot.sugar.AR     = MdlPS.AR{1};
pproc_boot.sugar.Sigma  = sqrt(MdlPS.Variance);


% Build outside option short series (this will be used to match actual decision):
% Select prices for years used in the estimation:
prices_used = [];
for y = [2003:2013] % Those are all available years of canasat data. opt.Years will be used later for selection in estimation
    prices_used = [prices_used; double(prices(prices.year == y, {opt.CornPrice opt.SoyPrice opt.SugarPrice}))];
end

% Get results from outside option model estimation:
if opt.WhichOutside == 1
    [OutDS, OutOptionGrid] = OutsideOptionModel(DS,opt);
elseif opt.WhichOutside == 2
    [OutDS, OutOptionGrid] = OutsideOptionModelSimple(DS,opt);
elseif opt.WhichOutside == 3
    [OutDS, OutOptionGrid] = OutsideOptionModelOld(DS,opt);
end

OutIndices = unique(OutDS.Index)';

ExoTimeStates        = nan(size(DS,1),size(prices_used,1));
RegionsForExoStates  = nan(size(DS,1),1);

for r = OutIndices
    
    if r ~= 0
    
    bootout(1) =  pproc.outside{r}.Cte/(1-pproc.outside{r}.AR);   
    for t=2:300
        bootout(t) = pproc.outside{r}.Cte + pproc.outside{r}.AR*bootout(t-1) + randn(1,1)*pproc.outside{r}.Sigma;
    end
    bootout = bootout(end-T+1:end);
    
    
    % Create OuReturn variable to be matched to data:
    OutReturn(:,r) = OutOptionGrid.Base(r) + OutOptionGrid.Partialcorn(r)*(prices_used(:,1) - OutOptionGrid.PriceAverage(1)) ...
        + OutOptionGrid.Partialsoy(r)*(prices_used(:,2) - OutOptionGrid.PriceAverage(2));
    OutReturn(:,r) = OutReturn(:,r)/1000;
    
%    disp([' Estimate process for class ' num2str(r) ' return: '])
    % Estimate AR(1) for outside option return:
    MdlOR = arima(1,0,0);
    MdlOR = estimate(MdlOR,bootout','Display','off');
    
    % Price process:
    pproc_boot.outside{r}.OutIndex = r;
    pproc_boot.outside{r}.Cte      = MdlOR.Constant; % This is the the constant.
    pproc_boot.outside{r}.AR       = MdlOR.AR{1};
    pproc_boot.outside{r}.Sigma    = sqrt(MdlOR.Variance);
    
    % Discretize the price process by region:
    
    Sigma = [pproc_boot.sugar.Sigma^2,  0; 0, pproc_boot.outside{r}.Sigma^2];
    F     = [pproc_boot.sugar.AR, 0; 0, pproc_boot.outside{r}.AR];
    F0    = [pproc_boot.sugar.Cte; pproc_boot.outside{r}.Cte];
    
    % Discretize process:
    [pproc_boot.P{r}, pproc_boot.y{r}, pproc_boot.x{r}, MVndx] = VarTauchen(opt.pproc.Ns,F0,F,opt.pproc.bandwidth,Sigma);
    
    % Now we need to classify prices in the data:
    SugOut = [prices_used(:,3), OutReturn(:,r)];
    
    Ny = size(pproc.y{r},1);
    Nd = size(SugOut,1);
    
    K1 = kron(SugOut(:,1),ones(1,Ny)) - kron(pproc.y{r}(:,1)',ones(Nd,1));
    K1 = abs(K1);
    
    K2 = kron(SugOut(:,2),ones(1,Ny)) - kron(pproc.y{r}(:,2)',ones(Nd,1));
    K2 = abs(K2);
    
    
    [~,exo_idx(:,r)] = min(K1+K2,[],2); 
    clear K1 K2
    
    NumInOut = sum(OutDS.Index == r);
    ExoTimeStates(OutDS.Index == r,:) = repmat(exo_idx(:,r)', [NumInOut,1]);
    RegionsForExoStates(OutDS.Index == r) = repmat(r, [NumInOut,1]);
    
    end
    
    % Table with exo time states:
    
    
end

end




